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New partial differential equations for the Wiener-Hermite expansions of the Langevin (stochastic) 
transitions are formulated. They are solved recursively in full order series solutions with respect 
to \ft. A sort of 'gauge' degrees of freedom (arbitrariness) involved in the solutions are analyzed 
and clarified. The Wiener-Hermite expansions play important roles as basic elements of numerical 
simulations of the Langevin equations. Specific solutions within some orders are presented as ex- 
amples. These expansions giving integral representations of the Fokker-Planck evolution kernels, 
similar formulations are possible for the imaginary time Hamiltonian evolution kernels as well. 

j< ". PACS 02.50.-r, 02.60.-x, 05.40.+j 

It is a long time since an idea that one expands stochastic variables of interest in terms of some normal ones, say, 
whose distributions are Gaussians, was emphasized by Wiener as an important concept when discussing stochastic 
processes, e.g., diffusions. Diffusion processes with drift are modeled in the Langevin equations in use of idealized 
white noise terms as, Xi = Ui{x) + rji(t) (i = 1, 2, • • • , d), with (f]i(t)Vj(t')) — 25y£(i — t'), and (r]i(t)) = 0, where 
Xi denotes a position of a particle in d-dimensional flat space. These stochastic differential equations define a 
I ■ transition, Xi(t) = Xi(t) — xoi for a finite time lapse t with an initial condition, 2^(0) = xqi- It is a natural question to 
ask if one can obtain a suitable expansion of the Langevin transition JQ (t) in terms of Gaussian random variates so 
that the expansion should have the same transition density as Xi(t), i.e., a weak solution. One calls this the Wiener- 
Hermite expansion The Langevin equations play important roles in wide range, from theoretical issues in 
physics such as stochastic quantization || together with its numerical applications in lattice field theories || , to many 
applications in various fields of science and technology . These expansions have proper significance in numerical 
simulations of the Langevin equations since a finite number of Gaussian variates replace white noise, naively, a set of 
infinite number of Gaussian variates. Since explicitly solvable Langevin equations |^,^| are restricted, approximations 
of solutions become important in practical applications. One can see in an extensive review |?J that quite a few 
approximation schemes have been examined in a long history of numerical applications of the Langevin equations, 
and yet the problem of systematic analysis of the approximation schemes is open. A reasonablly systematic analysis 
toward higher order approximation schemes (weak Taylor approximation scheme was reported lately within any 
■ fixed finite order flo|| . It is, however, welcome to obtain the Wiener-Hermite expansions of the Langevin transitions 
Xi(t) in full order series solutions with respect to \ft. Then these solutions will serve as basic elements for development 
<D ■ of new schemes as Runge-Kutta like ones. 

In this paper, we report a new successful formulation to obtain the Wiener-Hermite expansions in use of partial 
differential equations of Xi (t) with respect to \ft and d independent Gaussian random variates £j with variance 2 and 
mean 0. They are solved recursively in full order series solutions with respect to \fi. A sort of 'gauge' degrees of 
freedom (arbitrariness) can be involved in the solutions, and they are analyzed and clarified below. 

Formulation and solutions in the flat space. One can derive from the Langevin equations (see, e.g., |^-|^]) the 
corresponding Fokker-Planck equations dtP{t, x; 0, xq) = KP(t, x; 0, xo), for the transition density P, where K — 
di(di — Ui(x)), with an initial condition P(0, x; 0, Xo) = 6(x — xq). Now the problem to solve is to find a vector 
function Xi(t) whose variables are t and £i, with parameters, values of u^s and their higher derivatives evaluated at 
xq, satisfying 



> 



X 



e 



5(x-x ) = (5{x-x -X(t))) , (1) 



where (• • •) is averaging over the Gaussian distribution of This equation (1) can be rewritten as a sort of intergal 
representation of the Fokker-Planck evolution kernel as follows, 
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>K 'x ) = (x\(e- x -°)\x ) , (2) 



where Dirac's notation is used, and it is to be noted that K — di(di — Ui(x)) and that the evolution operator on 
the l.h.s. is arranged in normal order (all <9's in left of all x's). By differentiating (2) with respect to t, to obtain 
((A' + dtX ■ d x )e~ x ' dx )(x\xo) — , and then replacing the argument x of Ui in K with xq + X, and dropping off 
(a;|a;o), one derives an operator equation as, 

{{d 2 x -u-d x + d t x- d x ) e - x - 9 *) = , (3) 
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where one may consider either that only differential operators d x are involved with all the other parameters or that all 
operators are normal-ordered. Putting X{ = s^, where s = \ft, one assumes that the unknown ^ is to be expanded 
as £■ = £j+higher orders w.r.t. s, as a standard choice 0,|ffj- It implies that £' and £ can be expanded perturbatively 
with each other, and that the unknown is the coordinate transformation between them. Replacing d x in (3) with 
— (l/s)D', where and hereafter D' and D denote d^> and d^, respectively, and noting (...) = const x J d£ e~^ 2 
and making use of the partial integration formula and of change of integration variables to £' to have a form of the 
Fourier transformation, one obtains after some calculation on (3) the following partial differential equation as, 



(^)d^W - 5» - - *<«" + *>}] = ' (4 » 



where matrix notation is used, e.g., (D^')ij — d^'j, and summation over repeated indices should be understood. If 
one puts ^ = £i + syi , then one can derive from (4) the following partial differential equation as, 



(wL{ (1 + Sds + M)ijVj ~ 8 & ~ 2Dj) ( l + sDy {Dv)2 ) l3 ~ 2Ul{X0 + = ° ' (5) 



where an operator is defined as M^- = — 2Dj). Since we assume series solutions of yi with respect to 

s, and assume that each order is given by a polynomial of £ i; we define a space F of vector functions as F = 
{fi = pi e~i^ | pi is a polynomial of £. }. A basis of this space can be spaned by those functions f^s as, fi = 
Ci^i l ^ 2 ^^ n D il D i2 ■ ■ ■ D in e~i^ , where the polynomial part is essentially the generalized Hermite polynomial. It is 
to be noted that any coefficient tensor C^i^,.^ can be decomposed uniquely into totally symmetric component 
with respect to all indices including j and the residual component. The decomposition of coefficient tensors defines 
decomposition of F. We denote a space of the totally symmetric component as Fg, and a space of the residual 
component as Fp^s, and F = Fs + -Fpjvs- Here we name fi € F$ as a symmetric vector function, and fi € Fp^s 
as a pure nonsymmetric vector function. Then for fi G F, Difi = implies fi is a pure nonsymmetric vector 
function. A space P of polynomial vector function pi can be defined as P = {pi \ pi — e*^ fi, fi S F}, and 
P = Ps + Ppns with similarly defined Ps and Ppns, where Ps is a space of symmetric vector polynomials and Ppns 
pure nonsymmetirc vector polynomials. One sees for pi € Ppns, that MijPj = 0, and can show for pi € Ps, in a form 
as k = e^ 2 Ci 7 i lt i 2t . t i n D il D i2 ■ ■ ■ D in e~^ 2 , that M^pj = (n+ fli2| . Thus it is clear how P splits into eigenspaces 
of M. Now one can derive the following equation from (5) in use of an arbitrary € Ppns, a s, 

{(1 + sd s + M)(y -u- a)}. = S (& - 2D,-) ( T -I^-(A,) 2 )^ + 2 (e^ +s2 ^ - lju, + s((l + sd^Dy)^ , (6) 

where the Taylor expansion of Ui(xo + s£') at xo should be understood. Since eigenvalues of the operator 1 + sd s + M 
are positive in view of our assumption for solutions yi , one easily obtains the following recursion formula for a scries 
solution yi = Y^=o s ™2/i(n), for n > as 



WW = «.(»> + TT j TJJ {& - MtTs15y (DY ^ v + - s (^ +s2Yyd * - 1H + (d + WAWDY)*}^ 

(7) 



where = X)fc=o sk<1 i(k), an d ^ = X)fe=o sk Vi(k), an d {• ■ •}(»«) denotes m-th order terms of {• • •}, and for n = 0, 
2/j(o) — u i + a i(o) ■ The partial differential equation (4) and the recursion formula (7) for its series solution are main 
results of the present work. One notes that arbitrariness of the solution, a sort of 'gauge' degrees of freedom, comes 
from ai G Ppns, and it is interesting to note that e^e - ** appears as a transverse vector function in a space of 
d-dimensional Gaussian random variates. Actual calculation of 1/(1 + n + M) •{••■} in (7) is easy in principle, and it 
is readily done if decomposition of {• • •} in eigenfunctions of M is obtained. Let pi be a given fc-th order polynomial 
of £. Then pi is written as a linear combination of k + 2 unknown eigenfunctions. By multiple operations of M on 
the equation, one obtains sufficient number of independent equations to solve the unknowns. 

Example ofO(s e ) expansion of Xi in the flat space. Specific example of the solution Xi given by (7) can be chosen 
arbitrarily to some extent as stated above. We give an example in simplicity principle as follows, 

s 3 s 4 s 5 

Xi(s, f, x ) = s£i + s 2 Ui + —£,jdjUi + —^j^i 'ljih u, + 3ujdjUi + d 2 Ui) + — (S^jUkdjdk^ + A^djUkdkUi + i ; <V 4 // ,<},,//, 

s 6 

+ijikS, m d ] dkd m u t + 2£ j djd 2 u i ) + — {2djUkdjd k u t + UjUkdjdkU, + d UjdjUi + UjdjUkdkUt + 2u 1 djd 1 u i + d 2 d 2 u l ) , (8) 
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where, for brevity's sake, those eigenfunctions with eigenvalues other than M = 1 are dropped off in 0(s 6 ), which 
only affect higher orders in calculation of expectation values of any quantities. In the process of the above derivation, 
dj is so chosen as to cancel out djUj from the expression. The result (8) coincides with the former result obtained in 
a different approach jl9], except for allowable deviation of M = 4 eigenfunction in 0(s 5 ) terms for the same reason 
as in 0(s 6 ). We give a brief interpretation on that arbitrariness encountered here other than those in the full order 
expansion of Xj, which is characteristic to the finite order approximation. Suppose that one obtains the full order 
series solution Xj and one is to use it in the transition density (1) within some finite order, say 0(s 6 ). It is needless to 
say that one only has to retain terms in X up to 0(s 6 ), and it is only required to give correct expectation values of all 
moments of Xj within 0(s 6 ). Let us consider, for example, roles of 0(s 5 ) terms in calculation of ((Xj) m ), and we note 
that only relevant cases are when m = 1, 2, since the leading order of Xj is O(s) and that (0(s 5 ) terms) = if m = 1. 
When m = 2, only (s£ • (0(s 5 )terms)) is important. Key observation is that ((£) fe • (M = n)eigen function) = if 
k < n — 1. Thus (M = n)eigen function with n = 4, 6, • • ■ may be subtracted in 0(s 5 ) terms in Xj within the 0(s e ) 
approximation scheme. In the previous work |fjof , actual performance of the present 0(s 6 ) algorithm in numerical 
simulation of Langevin equation was tested in application to a single degree of freedom U(l) statistical model, where 
a significant improvement over 0(s 4 ) algorithm was observed in data in a theoretically expected manner. 

Formulation and solutions in the Lie group space. The present method for the Wiener-Hermite expansion can 
be extended to the Langevin transition in the Lie group space. We give a brief intepretation on the formulation 
and its solutions. Specifically we assume the Lie group is SU(N). The Langevin equation in the case Jf||,[l4| reads 
as g~ 1 -§£g = K(ui + T]i) where g S SU(N) and Aj's (i — 1, ■ • • , N 2 — 1) form antihermitean basis of Lie algebra, 
[Aj,Aj] = CijkXk with normalization tr(X^Xj) = Sij, and iji's are white noises as given above. The corresponding 
Fokker-Planck equation reads as dt4> = K<p, where 4> is normalized as J dfj, cj) = 1 with the Haar measure d/J,, and 
K = Vj(Vj — Ui), and Vj's are right differentiations as defined in e e V f{g) = f(ge eX ) for any functions / on SU(N), 
and it holds [Vj,Vj] = c^feVfc. Definition of the transition Xj is given as e tK S mv (g — go) = (e~ x ' v S mv (g — go)), 
where S mv is an invariant S-function, and it holds S inv (ge~ x — go) = S mv (g — goe x x ). Line of reasoning to reach 
a partial differential equation for Xj is parallel to that of (4) except for non-commutativity of differentiation Vj's. 
We make use of left differentiation T>i with respect to e x ' v such that X^e 1 = V^e 3 ^, where — TVs satisfy the same 
commutation relations as Aj's and can be written as X>j = [h{x ■ C))ijd x . with a function ft, of a single variable x, 
h(x) = (—x)/(e~ x — 1), and with matrices {Ci)jk = Cikj, [i = L ■ ■ ■ , N 2 — I). The corresponding equation to (4) 
reads as follows, 



D k 



where hij = (h(—X ■ C))ij and Xj = s££. To (5) corresponds the following equation as, 



-((h'L.h-L^l-sDy)) } 

S \ / ij J 



o 



(10) 



where Lj = ^ — 2Dj and h — I + sw are operators, and w = w(sy,^, s). Now it is derived in the same way as before 
in use of an arbitrary Oj € Ppns that 



{(1 + sd s +M){y-u-a)} i = sftfLjh 1 + \ D (By) 2 ) ^ + 2{tie^ +s2 ^ - l) lk u k + (rfLjW^ + w%i 

-si^utLj + LjW + sw t L j w)(DY)Y + s((l + sd s )a j )(Dy) :ji . 
The recursion formula for a series solution yi = X)5^Lo sn Vi(n)i f° r Ti > reads as 

Vi(n) = Oi(n) + 1 ~ ^ l DY (DY) 2 )^ + 2(if t e^+ sV > v - l) lfcUfc + ( J ff t £ i W)j j 

" + LjW + sW^W^DY)) + s((l + sdM^DY)^} , 

\ / y J (n) 



(11) 



(12) 



where ij = X)fe=o sk Di{k), — Sfc=o sk(1 i(k), W = w(sY,£, s), and = 1 + sVK on the right hand side, and for 
n = 0, 2/j( ) = ^ + a i( o). 
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Example of 0(s 4 ) expansion of Xj in the Lie group space. Specific example of the solution Xi due to (12) can be 
given as follows, 

s 3 / N \ s 4 

Xi(s, £, So) = s& + s 2 u 4 + — ^-VjUj - — CiJ + y (3V 2 u, + 3 Mj ■ + iVui) , (13) 

where for brevity's sake, only M = 1 eigenfunction is retained in 0(s 4 ) for the reason characteristic to finite order 
approximation as stated above, and u$ = V^S 1 is used. The result (13) shows that (12) reproduces correctly the one 



previously obtained 13 



Discussions and outlooks Some remarks are in order with respect to the Wiener- Hermite expansion of the Langevin 
transition Xi in the present formulation. First of all, in view of (1) and (2), it is obviously useful in numerical 
simulation of distribution (x\e TK \xo) with a fixed xq, and T = nt, where a choice of n sets of d independent Gaussian 
random varietes, i.e., those in nd dimensions, determines a discrete path, and an 0(t k+1 ) approximation of Xj yields 
0{l/n k ) one for the distribution for the time T Jl^Jt]]. Another point to remark is that one obtaines a distribution 
function (x\e tK \xo) of x by finding an inverse function £ = of £'(£) and by substituting £' = (x — x Q )/s to 

const ■ det(^r)e~3^ . A couterpart in the Lie group space of this observation holds as well, if one defines properly 
a distribution function t/j in noncompact space of the canonical coordinate of the first kind, X, with g = goe x x so 
that J* dX ip(X)f(g) = j d[i(g)(g\e tK \go)f(g), and then it is given as ip — const ■ det(J^-)e~4? \z'=x/s- In view of the 
defining equation (2) of the Wiener- Hermite expansion Xi, one can see obvious extensions of the present formulation 
to cases where on the left hand side of (2) the Fokker-Planck evolution operators are arranged in other orderings, e.g., 
the Weyl ordering, instead of the normal one. Then one can have the Xi depending on both parameters xq and x, 
which is so called implicit algorithms in j?J . Another extenstion of the present formulation is toward the case where 
the imaginary time Hamiltonian evolution operators e~ tH are treated in various operator orderings. Then we need 
an extra scalar function in the exponent of the right hand side of (2), since the distribution is not normalized any 
more. This extension will also give means of numerical simulations of the ground state in such systems of interest 
provided shrinking of the distribution is properly taken care of. Finally we comment that a ^-linear solution is 
correctly reproduced from (7) with = for the Ornstein-Uhlenbeck process ||, the case itj being linear functions. 
We conclude that main results of the present work are the partial differential equations (4) and (9) together with 
their solutions (7) and (12), respectively, and that although their mathematical rigorousness on a level of such 
as convergence of the series solutions is beyond scope of the present analyses, we believe those will play important 
roles as a basis for systematic development of weak approximation schemes including Runge Kutta like ones. 

The author is grateful to S. Furui for constant colaboration and discussions on numerical simulations from which 
motivation of the present analyses came up, and he is also grateful to Y. Oyanagi for bringing to his attention the 
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